3d ball skinning using partial differential equations for generation of smooth tubular surfaces

ABSTRACT

A method of computing a continuous interpolation of a discrete set of three-dimensional (3D) balls, including generating an initial skin, wherein the initial skin is a surface comprised of splines and wherein the splines touch each ball along a circle that is tangent to the ball, solving a first differential equation to minimize the initial skin&#39;s surface area or solving a second differential equation to minimize a squared mean curvature of the initial skin&#39;s surface, wherein the result of solving the first or second differential equations is an updated skin; and repeating the steps of solving the first or second differential equations for the updated skin, and then, repeating the steps of solving the first or second differential equations for each subsequently updated skin until a desired skin is realized.

CROSS-REFERENCE TO RELATED APPLICATION

This application claims the benefit of U.S. Provisional Application No.61/091,494, filed Aug. 25, 2008, the disclosure of which is incorporatedby reference herein in its entirety.

BACKGROUND OF THE INVENTION

1. Technical Field

The present invention relates to three-dimensional (3D) ball skinning.

2. Discussion of the Related Art

The geometric problem of ball skinning is the computation of acontinuous interpolation of a discrete set of balls, for example. Thisproblem arises in numerous applications, including character skinning,molecular surface model generation, and modeling of tubular structures,etc.

In computer animation, often an articulated object or character isconstructed using a layered representation consisting of a skeletalstructure and a corresponding geometric skin, see e.g., K. Singh, E.Kokkevis, Skinning Characters using Surface Oriented Free-FormDeformations, in: Graphics Interface, 2000, pp. 35-42. Here, theskeleton has fewer degrees of freedom and is simpler to adjust by ananimator. Given a new skeletal pose, the skinning algorithm isresponsible for deforming the geometric skin to respond to the motion ofthe underlying skeleton. The skinning problem is also a special case ofthe problem of computing the envelopes of families of quadrics, whichhave been investigated by Peternell in, M. Peternell, RationalParameterizations for Envelopes of Quadric Families, Ph.D. thesis,University of Technology, Vienna, Austria (1997), via the use ofcyclographic maps. Rossignac and Schaefer in, J. Rossignac, S. Schaefer,J-splines, Computer Aided Design 40 (10-11) (2008) 1024-132, presentedJ-splines, which produce smooth curves from a set of ordered pointsusing a subdivision framework.

The problem of ball skinning appears frequently in the context ofcomputational chemistry and molecular biology, when generating surfacemeshes for molecular models, see e.g., H. Cheng, X. Shi, Quality MeshGeneration for Molecular Skin Surfaces Using Restricted Union of Balls,in: IEEE Visualization, 2005, H. Edelsbrunner, Deformable smooth surfacedesign, Discrete and Computational Geometry 21 (1) (1999) 87-115 and N.Kruithov, G. Vegter, Envelope Surfaces, in: Annual Symposium onComputational Geometry, 2006, pp. 411-420. Several algorithms exist toskin a molecular model to produce a C¹ continuous surface that istangent smooth and has high mesh quality. These methods are typicallyeither based on Delaunay triangulation like Cheng's or by finding theisosurface of an implicit function like Kruithov's. The work of Kruithovet al. derives a special subset of skins that is piece-wise quadratic.When dealing with a continuous family of balls, the skin may be computedas the envelope of the infinite union of the circles of intersection oftwo consecutive balls of infinitely close center. While the surfacesgenerated by these methods are tangent to the balls and have smoothnessat the point of tangency, none of these existing methods provides anoptimally smooth skin.

SUMMARY OF THE INVENTION

In an exemplary embodiment of the present invention, a method ofcomputing a continuous interpolation of a discrete set ofthree-dimensional (3D) balls, comprises: generating an initial skin,wherein the initial skin is a surface comprised of splines and whereinthe splines touch each ball along a circle that is tangent to the ball;solving a first differential equation to minimize the initial skin'ssurface area or solving a second differential equation to minimize asquared mean curvature of the initial skin's surface, wherein the resultof solving the first or second differential equations is an updatedskin; and repeating the steps of solving the first or seconddifferential equations for the updated skin, and then, repeating thesteps of solving the first or second differential equations for eachsubsequently updated skin until a desired skin is realized.

A circle on a ball B^(i) is in a plane with normal N^(i)=[cos θ^(i) sinφ^(i), sin φ^(i) sin φ^(i), cos φ^(i)]^(T) that passes through the ballcenter c^(i) and intersects the ball B^(i), q(u) is a point on thecircle and is defined by q(u)={circumflex over (x)} cos u+ŷ sin u, and apoint p^(i)(u) on a circle through which a spline passes is defined byp^(i)(u)=c^(i)+r^(i)R^(i)q(u), where r^(i) is a radius of an ith ball,and R^(i) is a rotation matrix specified by N^(i).

The rotation matrix R^(i) is

$R^{i} = {{R\left( {\theta^{i},\varphi^{i}} \right)} = {\begin{bmatrix}{\cos \; \theta^{i}\cos \; \varphi^{i}} & {{- \sin}\; \theta^{i}} & {\cos \; \theta^{i}\sin \; \varphi^{i}} \\{\sin \; \theta^{i}\cos \; \varphi^{i}} & {\cos \; \theta^{i}} & {\sin \; \theta^{i}\sin \; \varphi^{i}} \\{{- \sin}\; \varphi^{i}} & 0 & {\cos \; \varphi^{i}}\end{bmatrix}.}}$

The initial skin is represented by S(u,v), where S(u,v) is a collectionof continuous segments S^(i)(u,v), where an S^(i)(u,v) is specified byA^(i)(u),B^(i)(u),C^(i)(u),D^(i)(u), whereA^(i)(u)=−2p^(i+1)(u)+2p^(i)(u)+t^(i)N^(i)+t^(i+1)N^(i+1),B^(i)(u)=3p^(i+1)(u)−3p^(i)(u)−2t^(i)N^(i)−t^(i+1)N^(i+1),C^(i)(u)=t^(i)N^(i), and D^(i)(u)=p^(i)(u)

The first differential equation is represented by

${\frac{\partial J^{a}}{\partial w^{k}} = {{{\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2F^{i}\frac{\partial F^{i}}{\partial w^{k}}}}{\left\lbrack {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{v}}}}} + {\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2F^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}}}{\left\lbrack {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{{v}.\mspace{79mu} J^{a}}}}}}} = {\int{\int{\sqrt{{EG} - F^{2}}{u}{v}}}}}},{w^{k} \in {\left\lbrack {\theta^{i},\varphi^{i}} \right\rbrack {\forall i}}},\mspace{79mu} {E = {S_{u} \cdot S_{u}}},{F = {S_{u} \cdot S_{v}}},\mspace{79mu} {and}$     G = S_(v) ⋅ S_(v).

The second differential equation is represented by

$\frac{\partial J^{c}}{\partial w^{k}} = {{\int{\int{2\; {H^{i}\begin{bmatrix}{\frac{\begin{matrix}{{\frac{\partial e^{i}}{\partial w^{k}}G^{i}} + {e^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial f^{i}}{\partial w^{k}}F^{i}} -} \\{{2\; f^{i}\frac{\partial F^{i}}{\partial w^{k}}} + {\frac{\partial g^{i}}{\partial w^{k}}E^{i}} + {g^{i}\frac{\partial E^{i}}{\partial w^{k}}}}\end{matrix}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)} -} \\\frac{\begin{matrix}{\left( {{e^{i}G^{i}} - {2\; f^{i}F^{i}} + {g^{i}E^{i}}} \right) \cdot} \\\left( {{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial F^{i}}{\partial w^{k}}F^{i}}} \right)\end{matrix}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)^{2}}\end{bmatrix}}{u}{v}}}} + {\int{\int{2\; {H^{i - 1}\begin{bmatrix}{\frac{\begin{matrix}{{\frac{\partial\text{?}^{i - 1}}{\partial w^{k}}G^{i - 1}} + {e^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial f^{i - 1}}{\partial w^{k}}F^{i - 1}} -} \\{{2\; f^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}} + {\frac{\partial g^{i - 1}}{\partial w^{k}}E^{i - 1}} + {g^{i - 1}\frac{\partial E^{i - 1}}{\partial w^{k}}}}\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)} -} \\\frac{\begin{matrix}{\left( {{e^{i - 1}G^{i - 1}} - {2\; f^{i - 1}F^{i - 1}} + {g^{i - 1}E^{i - 1}}} \right) \cdot} \\\left( {{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial F^{i - 1}}{\partial w^{k}}F^{i - 1}}} \right)\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)^{2}}\end{bmatrix}}{u}{{v}.\text{?}}\text{indicates text missing or illegible when filed}}}}}$

J^(c)=∫∫H²dudv, w^(k)ε[θ^(i),φ^(i)]∀i, E=S_(u)·S_(u), F=S_(u)·S_(v),G=S_(v)·S_(v), e=M·S_(uu), f=M·S_(uv), and G=M·S_(vv), where M is thesurface normal.

The method further comprises combining the first and second differentialequations to generate the updated skin, wherein the first and seconddifferential equations are combined by the following equation

${\frac{\partial J}{\partial w^{k}} = {{\left( {1 - k} \right)\frac{\partial J^{a}}{\partial w^{k}}} + \frac{\partial J^{c}}{\partial w^{k}}}},$

where

$\frac{\partial J^{a}}{\partial w^{k}}$

is the first differential equation and

$\frac{\partial J^{c}}{\partial w^{k}}$

is the second differential equation.

The first and second differential equations are used in a gradientdescent procedure to find the desired skin.

The gradient descent procedure manipulates parameters w^(i)=[θ^(i),φ^(i)]^(T) of each ball i, where w^(i)(n+1)=w^(i)(n)−Δt∇J_(w) _(i)_((n).)∀i to find the desired skin.

In an exemplary embodiment of the present invention, a method ofmodeling a tubular structure, comprises: imaging a tubular structure;placing a plurality of balls in the tubular structure; and finding askin that smoothly interpolates the balls, wherein finding the skincomprises: generating an initial skin, wherein the initial skin is asurface comprised of splines and wherein the splines touch each ballalong a circle that is tangent to the ball; solving a first differentialequation to minimize the initial skin's surface area or solving a seconddifferential equation to minimize a squared mean curvature of theinitial skin's surface, wherein the result of solving the first orsecond differential equations is an updated skin; and repeating thesteps of solving the first or second differential equations for theupdated skin, and then, repeating the steps of solving the first orsecond differential equations for each subsequently updated skin untilthe skin that smoothly interpolates the balls is realized, wherein theskin that smoothly interpolates the balls is a model of the tubularstructure.

The tubular structure is an anatomical structure.

The tubular structure is imaged by a scanner.

The balls are ordered.

In an exemplary embodiment of the present invention, a system ofcomputing a continuous interpolation of a discrete set of 3D balls,comprises: a memory device for storing a program; and a processor incommunication with the memory device, the processor operative with theprogram to perform a method, the method comprising: generating aninitial skin, wherein the initial skin is a surface comprised of splinesand wherein the splines touch each ball along a circle that is tangentto the ball; solving a first differential equation to minimize theinitial skin's surface area or solving a second differential equation tominimize a squared mean curvature of the initial skin's surface, whereinthe result of solving the first or second differential equations is anupdated skin; and repeating the steps of solving the first or seconddifferential equations for the updated skin, and then, repeating thesteps of solving the first or second differential equations for eachsubsequently updated skin until a desired skin is realized.

The foregoing features are of representative embodiments and arepresented to assist in understanding the invention. It should beunderstood that they are not intended to be considered limitations onthe invention as defined by the claims, or limitations on equivalents tothe claims. Therefore, this summary of features should not be considereddispositive in determining equivalents. Additional features of theinvention will become apparent in the following description, from thedrawings and from the claims.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is ball skinning according to an exemplary embodiment of thepresent invention;

FIG. 2 is a representation of a circle of contact according to anexemplary embodiment of the present invention;

FIG. 3 shows segments and a spline according to an exemplary embodimentof the present invention;

FIG. 4 is ball skinning according to an exemplary embodiment of thepresent invention;

FIG. 5 is ball skinning according to an exemplary embodiment of thepresent invention;

FIG. 6 is ball skinning according to an exemplary embodiment of thepresent invention;

FIG. 7 is ball skinning according to an exemplary embodiment of thepresent invention;

FIG. 8 is a convergence plot for the exemplary embodiment shown in FIG.7;

FIG. 9 is a comparison between ball skinning according to an exemplaryembodiment of the present invention and J-splines;

FIGS. 10A and 10B illustrate convergence to desired and undesired minimabased on poor and severely poor initializations in ball skinningaccording to an exemplary embodiment of the present invention;

FIG. 11 illustrates skin passing through a ball in ball skinningaccording to an exemplary embodiment of the present invention;

FIG. 12 is a flowchart illustrating a method of ball skinning accordingto an exemplary embodiment of the present invention; and

FIG. 13 is a block diagram of a system in which exemplary embodiments ofthe present invention may be implemented.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

In an exemplary embodiment of the present invention, we are interestedin modeling the geometry of a blood vessel that has been identifiedusing Pearling as described by B. Whited, J. Rossignac, G. Slabaugh, T.Fang, G. Unal, Pearling: 3D Interactive Extraction of Tubular Structuresfrom Volumetric Images, in: Interaction in Medical Image Analysis andVisualization, held in conjunction with MICCAI, 2007. Pearling is a ballpacking algorithm that places numerous balls of different radii so thatthey fit snugly inside an imaged vessel. Given these balls, we wouldlike to find a C¹ skin that smoothly interpolates the balls. Thissurface can then be used for visualization of the blood vessel,simulation of blood flows using computational fluid dynamics, as well asmeasurements such as volume or surface area. We note that the problem oftwo-dimensional (2D) ball skinning was addressed in our paper, G.Slabaugh, G. Unal, T. Fang, J. Rossignac, B. Whited, VariationalSkinning of an Ordered Set of Discrete 2D Balls, in: Geometric Modelingand Processing, 2008 and U.S. Patent Application Publication No.2008/0109192, entitled “System and method for variational ball skinningfor geometric modeling of ordered balls”, the disclosures of which areincorporated by reference herein in their entirety. In the paper, forexample, we extended the methodology to skinning three-dimensional (3D)balls. In accordance with the present invention, the problem has asimilar conceptual formulation based on differential equations; however,the geometry is notably different: instead of minimizing the arc lengthand curvature of a curve, we minimize the surface area and meancurvature of a surface. Furthermore, the geometry of contact between theskin and the balls, as well as the differential geometry of the skin issignificantly different in the 3D case.

In accordance with the present invention, our approach produces asurface that minimizes a convex combination of surface area and squaredmean curvature terms. Such surfaces are loosely referred to as “minimalsurfaces” in the literature. It has been shown in the surface evolutionliterature that, for unconstrained surfaces, these terms result incurvature and Willmore flows, see e.g., J. A. Sethian, Level Set Methodsand Fast Marching Methods Evolving Interfaces in Computational Geometry,Fluid Mechanics, Computer Vision, and Materials Science, CambridgeUniversity Press, 1999 and K. E., R. Schatzle, Gradient Flow for theWillmore functional, Comm. Anal. Geom. 10 (5), respectively. By addingthe additional constraint that the skin must pass through a circle ofintersection with each ball, we significantly reduce the dimensionalityof the optimization.

In a method of ball skinning according to an exemplary embodiment of thepresent invention, we model the skin as a C¹ surface, which, byconstruction, must touch each ball along a circle that is tangent to theball. We then provide two novel derivations, one for deforming thisconstrained surface to minimize its surface area; and a secondderivation for minimizing its squared mean curvature. The result ofthese derivations is differential equations, which we then solve toupdate a given surface to its optimal position. We then showexperimental examples of how these differential equations are used toperform optimally smooth skinning of balls.

In accordance with an exemplary embodiment of the present invention,with regard to the representation of a skin's surface, we require theskin to pass through a circle on the ith ball B¹, as depicted in FIG. 2(a). This circle resides in a plane with normal N¹=[cos θ^(i) sin φ^(i),sin φ^(i), sin φ^(i) sin φ^(i), cos φ^(i)]^(T) that passes though theball center c^(i) and intersects B^(i).

Let q(u) be a point on the unit circle in the z=0 plane, as shown inFIG. 2( b). We can parameterize q(u) as

q(u)={circumflex over (x)} cos u+ŷ sin u,   (1)

where u is a parameterization angle. We can then express p^(i)(u), apoint on the circle of B^(i) as

p ^(i)(u)=c ^(i) +r ^(i) R ^(i) q(u),   (2)

where r^(i) is the radius of the ith ball, and R^(i) is a rotationmatrix specified by N^(i) as

$\begin{matrix}{R^{i} = {{R\left( {\theta^{i},\varphi^{i}} \right)} = \begin{bmatrix}{\cos \; \theta^{i}\cos \; \varphi^{i}} & {{- \sin}\; \theta^{i}} & {\cos \; \theta^{i}\sin \; \varphi^{i}} \\{\sin \; \theta^{i}\cos \; \varphi^{i}} & {\cos \; \theta^{i}} & {\sin \; \theta^{i}\sin \; \varphi^{i}} \\{{- \sin}\; \varphi^{i}} & 0 & {\cos \; \varphi^{i}}\end{bmatrix}}} & (3)\end{matrix}$

p^(i)(u) provides a way to parameterize the circle on each ball. Wewould like to form a parametric skin S(u,v) that satisfies severalgeometric criteria:

-   -   (1) The skin should be modeled by a circle that contacts each        ball.    -   (2) The skin should be tangent to each ball along the circle of        contact.    -   (3) The skin should optimize an energy functional composed of        surface area and mean curvature terms.

We compose the skin S(u,v) as a set of segments S^(i)(u,v) for i=1 . . .N−i, where N is the total number of balls, as depicted in FIG. 3. Toform the segment S^(i)(u,v), we would like to generate a spline-basedsurface that connects the circles on adjacent balls. Various splinerepresentations (such as Catmull-Rom, 4-point, etc.) are possible formodeling segments using a set of splines. Each spline starts at pointp^(i)(u) in direction N^(i), and ends at point p^(i+1)(u) in directionN^(i+1), with

S ^(i)(u,v)=A ^(i)(u)v ³ +B ^(i)(u)v ² +C ^(i)(u)v+D ^(i)(u),   (4)

since the four constraints require four degrees of freedom. For the ithsegment, A^(i)(u),B^(i)(u),C^(i)(u), and D^(i)(u) are coefficients, andv ε[0,1] is a parameterization variable.

Each segment S^(i)(u,v) of the skin is defined by the Hermiteinterpolation of the boundary conditions, specifically,

${{{S^{i}\left( {u,v} \right)}_{v = 0}} = {p^{i}(u)}},{{\frac{\partial{S^{i}\left( {u,v} \right)}}{\partial v}_{v = 0}} = {t^{i}N^{i}}},{{{S^{i}\left( {u,v} \right)}_{v = 1}} = {p^{i + 1}(u)}},{and}$${{\frac{\partial{S^{i}\left( {u,v} \right)}}{\partial v}_{v = 1}} = {t^{i + 1}N^{i + 1}}},$

where, for each i, t^(i) is a stiffness coefficient that controls theinfluence of the normal N^(i). Each t^(i) is fixed to be half thedistance between the next and previous ball centers (for the first andlast balls, it is the distance between the ball center and its neighborball center) for all examples in this disclosure. Such a stiffnessencourages smoothness of the connecting segments at a circle, and isbased on the central difference approximation of the first derivativecomputed using ball centers. We note however that straighter segmentscan be achieved by scaling the stiffness by a coefficient less than one.

With these constraints, and the derivative of the segment,

$\begin{matrix}{{\frac{\partial{S^{i}\left( {u,v} \right)}}{\partial v} = {{3\; {A^{i}(u)}v^{2}} + {2\; {B^{i}(u)}v} + {C^{i}(u)}}},} & (5)\end{matrix}$

we obtain a system of four equations for the four coefficients:D^(i)(u)=p^(i)(u),C^(i)(u)=t^(i)N^(i),A^(i)(u)+B^(i)(u)+C^(i)(u)+D^(i)(u)=p^(i+1)(u),and 3A^(i)(i)+2B^(i)(u)+C^(i)(u)=t^(i+1)N^(i+1), which is easily solved,yielding

A ^(i)(u)=−2p ^(i+1)(u)+2p ^(i)(u)+t ^(i) N ^(i) +t ^(i+1) N ^(i+1)

B ^(i)(u)=3p ^(i+1)(u)−3p ^(i)(u)−2t ^(i) N ^(i) −t ^(i+1) N ^(i+1)

C ^(i)(u)=t ^(i) N ^(i)

D ^(i)(u)=p ^(i)(u)   (6)

Thus, we have a way of describing the skin S(u,v) as a collection of N−1C¹ continuous segments S^(i)(u,v). Segment S^(i)(u,v) in turn isspecified by the coefficients A^(i)(u),B^(i)(u),C^(i)(u),D^(i)(u), andthese coefficients are functions of the circles p^(i)(u),p^(i+1)(u), andthe normals N^(i) and N^(i+1). Finally, the circles and normals are, inturn, functions of the angles φ^(i),φ^(i+1),φ^(i),φ^(i+1).

In accordance with an exemplary embodiment of the present invention,with regard to surface area minimization, we derive differentialequations to evolve the parameters of the skin to minimize the skin'ssurface area.

From differential geometry, it is well known that the surface area isgiven by

J ^(a)=∫∫√{square root over (EG−F²)}dudv,   (7)

where

E=S _(u) ·S _(u)   (8)

F=S _(u) ·S _(v)   (9)

G=S _(v) ·S _(v)   (10)

are coefficients of the first fundamental form.

Since S(u,v) is expressed as a sum of N−1 segments, this is equivalentto

$\begin{matrix}{J^{a} = {\sum\limits_{i = 1}^{N - 1}{\int{\int{\sqrt{{E^{i}G^{i}} - \left( F^{i} \right)^{2}}{u}{{v}.}}}}}} & (11)\end{matrix}$

We take the derivative of Equation 11 with respect to the parameterw^(k), where w^(k) ε[θ^(i),φ^(i)]∀i. This will give us a gradientdirection that we can use in a numerical gradient descent procedure tofind the angles that minimize the surface area of the skin. Since θ^(i)and φ^(i) affect only the ith and i−1st segments, we can replace thesummation with two surface integrals,

$\begin{matrix}{\frac{\partial J^{a}}{\partial w^{k}} = {{\frac{\partial}{\partial w^{k}}{\int{\int{\sqrt{{E^{i}G^{i}} - \left( F^{i} \right)^{2}}{u}{v}}}}} + {\frac{\partial}{\partial w^{k}}{\int{\int{\sqrt{{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}}{u}{{v}.}}}}}}} & (12)\end{matrix}$

Propagating the derivative through the integrals gives

$\begin{matrix}{\frac{\partial J^{a}}{\partial w^{k}} = {{\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\; F^{i}\frac{\partial F^{i}}{\partial w^{k}}}}{\left\lbrack {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{v}}}}} + {\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\; F^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}}}{\left\lbrack {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{{v}.}}}}}}} & (13)\end{matrix}$

Derivations for the derivatives of the coefficients of the firstfundamental form with respect to w^(k) are provided near the end of thisdisclosure.

In accordance with an exemplary embodiment of the present invention,with regard to curvature minimization, we derive differential equationsfor updating the skin to minimize its curvature. In 3D there are severalpotential curvatures one could employ, including mean and Gaussiancurvatures. In this disclosure, we focus on the mean curvature, which isclosely related to the first variation of surface area. The meancurvature is given by

$\begin{matrix}{{H = \frac{{eG} - {2\; {fF}} + {gE}}{2\left( {{EG} - F^{2}} \right)}},} & (14)\end{matrix}$

where E,F, and G are given in the previous section, and e,f, and g comefrom the fundamental form,

e=M·S _(uu)   (15)

f=M·S _(uv)   (16)

g=M·S _(vv)   (17)

where M is the surface normal (not to be confused with N, which is thenormal of the plane that intersects a ball).

Our energy is

J^(c)=∫∫H²dudv.   (18)

Since S(u,v) is expressed as a sum of N−1 segments, this is equivalentto

$\begin{matrix}{J^{c} = {\sum\limits_{i = 1}^{N - 1}{\int{\int{\left( H^{i} \right)^{2}{u}{{v}.}}}}}} & (19)\end{matrix}$

As before, we would like to take the derivative of Equation 19 withrespect to the parameter w^(k), where w^(k) ε[θ^(i),φ^(i)]∀i. This willgive us a gradient direction which we can use in a numerical gradientdescent procedure to find the angles that minimize the curvature of theskin. Since θ^(i) and φ^(i) affect only the ith and i−1st segments, wecan replace the summation with two surface integrals,

$\begin{matrix}{\frac{\partial J^{c}}{\partial w^{k}} = {{\frac{\partial\;}{\partial w^{k}}{\int{\int{\left\lbrack \frac{{e^{i}G^{i}} - {2\; f^{i}F^{i}} + {g^{i}E^{i}}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)} \right\rbrack^{2}{u}{v}}}}} + {\frac{\partial\;}{\partial w^{k}}{\int{\int{\left\lbrack \frac{{e^{i - 1}G^{i - 1}} - {2\; f^{i - 1}F^{i - 1}} + {g^{i - 1}E^{i - 1}}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)} \right\rbrack^{2}{u}{{v}.}}}}}}} & (20)\end{matrix}$

Propagating the derivative through the integrals gives

$\begin{matrix}{\frac{\partial J^{c}}{\partial w^{k}} = {{\int{\int{2\; {H^{i}\left\lbrack {\frac{\begin{matrix}{{\frac{\partial ^{}}{\partial w^{k}}G^{i}} + {e^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial f^{i}}{\partial w^{k}}F^{i}} -} \\{{2\; f^{i}\frac{\partial F^{i}}{\partial w^{k}}} + {\frac{\partial g^{i}}{\partial w^{k}}E^{i}} + {g^{i}\frac{\partial E^{i}}{\partial w^{k}}}}\end{matrix}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)} - \frac{\begin{matrix}{\left( {{e^{i}G^{i}} - {2\; f^{i}F^{i}} + {g^{i}E^{i}}} \right) \cdot} \\\left( {{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial F^{i}}{\partial w^{k}}F^{i}}} \right)\end{matrix}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)^{2}}} \right\rbrack}{u}{v}}}} + {\int{\int{2\; {H^{i - 1}\left\lbrack {\frac{\begin{matrix}{{\frac{\partial\text{?}^{i - 1}}{\partial w^{k}}G^{i - 1}} + {\text{?}^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial f^{i - 1}}{\partial w^{k}}F^{i - 1}} -} \\{{2\; f^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}} + {\frac{\partial g^{i - 1}}{\partial w^{k}}E^{i - 1}} + {g^{i - 1}\frac{\partial E^{i - 1}}{\partial w^{k}}}}\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)} - \frac{\begin{matrix}{\left( {{e^{i - 1}G^{i - 1}} - {2\; f^{i - 1}F^{i - 1}} + {g^{i - 1}E^{i - 1}}} \right) \cdot} \\\left( {{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial F^{i - 1}}{\partial w^{k}}F^{i - 1}}} \right)\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)^{2}}} \right\rbrack}{u}{{v}.\text{?}}\text{indicates text missing or illegible when filed}}}}}} & (21)\end{matrix}$

where the derivatives of the coefficients from the first fundamentalform (i.e., E, F, and G) are given near the end of this disclosure andthe derivatives of the coefficients from the second fundamental form(i.e., e, f and g ) are analytically derived near the end of thisdisclosure.

In implementing a method according to an exemplary embodiment of thepresent invention, we combine the energies J^(a) and J^(c) together, as

J=(1−k)J ^(a) +kJ ^(c),   (22)

where k is a constant used to weight the surface area minimizationrelative to the curvature minimization. Convex combinations of the twocan be selected using k ε[0,1]. Therefore, the combined energyminimization is given by

$\begin{matrix}{{\frac{\partial J}{\partial w^{k}} = {{\left( {1 - k} \right)\frac{\partial J^{a}}{\partial w^{k}}} + \frac{\partial J^{c}}{\partial w^{k}}}},} & (23)\end{matrix}$

where

$\frac{\partial J^{a}}{\partial w^{k}}$

is given in Equation 13 and

$\frac{\partial J^{c}}{\partial w^{k}}$

is provided in Equation 21. In all of the experiments in thisdisclosure, we fix k=0.9, to encourage smoother solutions.

These equations are a set of differential equations that can be used ina gradient descent procedure to optimize the skin by manipulating theparameters w^(i)=[θ^(i),φ^(i)]^(T) of each ball i. Let w^(i)(n) be theith ball's parameters at iteration n. We can then update the parametersby moving them in the negative gradient direction, i.e.,

w ^(i)(n+1)=w ^(i)(n)−Δt∇J _(w) _(i) _((n).∀) i,   (24)

where Δt is a time step.

The computational complexity of the algorithm depends on the number ofballs N and the number of points on the surface where the points andderivatives are evaluated. The number of points is given by LM, where Lis the number of sampling points on each spline, and M is the number ofsplines on a segment. For each iteration of the gradient descentprocedure, the computational complexity is O(NLM). The number ofiterations required depends on the time step Δt as well as how close theinitial skin is to the final solution.

A simple demonstration of ball skinning according to an exemplaryembodiment of the present invention is provided in FIG. 4. Here, fourballs of radius 3, 2, 2, and 3 units, respectively, were set in thexy-plane, at points c¹=[0,0,0]^(T), c²=[5,0,5]^(T), c³=[10,0,0]^(T) andc⁴=[17,0,−5]^(T). The initial parameters for this experiment were

${w^{1} = \left\lbrack {0,0} \right\rbrack^{T}},{w^{2} = \left\lbrack {0,\frac{\pi}{4}} \right\rbrack^{T}},{w^{3} = \left\lbrack {0,\frac{\pi}{2}} \right\rbrack^{T}},{{{and}\mspace{14mu} w^{4}} = \left\lbrack {0,\frac{\pi}{4}} \right\rbrack^{T}},$

respectively; the initial skin is shown in FIG. 4( a). The parameterswere iteratively updated using Equation 24, with L=50 and M=20 (thesevalues for L and M are used for all experiments in this disclosure). Anintermediate solution after 20 iterations is shown in FIG. 4( b); atthis stage, the skin is considerably smoother while satisfying theconstraints of the problem. We show the result after 40 iterations inFIG. 4( c), at which point the energy has reached a minimum and theparameters have converged. The solution (all 40 iterations) is computedin 4.3 seconds using C++ code on a machine with a 2.0 GHz processor. Werender the surface as a collection of splines, and additionally show thecircle of intersection on each ball. The energy of the surface, asmeasure using Equation 22, drops from 7.36×10⁹ in FIG. 4( a) to 670 inFIG. 4( c).

FIG. 5 is a slightly more complicated demonstration of ball skinningaccording to an exemplary embodiment of the present invention in whichsome balls overlap and others do not. The initial skin in shown in FIG.5( a), an intermediate result after 25 iterations in FIG. 5( b), and thefinal result upon convergence after 50 iterations in FIG. 5( c). Thesolutions (all 50 iterations) is computed in 10.4 seconds. The energy ofthe surface, as measured using Equation 22, drops from 1.86×10¹⁰ in FIG.5( a) to 1400 in FIG. 5( c).

FIG. 6 is a demonstration of ball skinning according to an exemplaryembodiment of the present invention for a symmetric configuration ofballs, but asymmetric initial conditions. The initial skin is shown inFIG. 6( a), an intermediate result after 30 iterations in FIG. 6( b),and the final result upon convergence after 60 iterations in FIG. 6( c).The solution (all 60 iterations) is computed in 8.3 seconds. The energyof the surface, as measured using Equation 22, drops from 2.79×10⁷ inFIG. 6( a) to 1287 in FIG. 6( c). Note that due to the symmetry of theballs, the skin itself is symmetric upon convergence.

More demonstrations of ball skinning according to an exemplaryembodiment of the present invention are provided in FIGS. 7 and 1.Initialization is shown in FIG. 7( a), an intermediate result after 25iterations in FIG. 7( b) and a final result upon convergence after 50iterations in FIG. 7( c). In FIG. 7, the balls are arranged on a sinewave and have a variable radius. In addition, some of the balls overlapwhile others do not. Convergence of the skinning algorithm, startingfrom a set of angles far from the optimal result, takes 11.4 seconds,and reduces the energy from 9.67×10⁵ to 11,471. In FIG. 1, the variableradius balls are arranged in a spiral (see FIG. 1( a)). The skin (seeFIG. 1( b)) is generated in 15.3 seconds.

In FIG. 8, we show a plot of the energy J of the surface vs. theiteration number. Note that initially, the energy is high and successiveiterations reduce the energy until convergence occurs around the 30^(th)iteration. Upon convergence, the energy oscillates around its minimalvalue. This fact can be exploited as an automatic convergence criterion.

We have implemented the J-splines technique of Rossignac and Schaefer towhich we compare a ball skinning method according to an exemplaryembodiment of the present invention. This approach outlines a generalsubdivision algorithm for producing smooth curves (J-splines) from a setof ordered points. Iterative applications of the subdivision algorithmyield a family of limit curves, one of which is a quintic b-spline (C⁴).This quintic b-spline will no longer interpolate the input points, butcan be “retrofitted” as Rossignac and Schaefer describes. Thisretrofitting process iteratively offsets the input control points untilthe final curve interpolates the input, thus resulting in aninterpolating C⁴ spline. As a comparison, we used this subdivisionapproach to subdivide the series of balls as a four-dimensional (4D)curve (x, y, z+radius). The cross sections of the skin are then computedas circles which lie on the surface of the convex hull of everyconsecutive pair of balls and are also orthogonal to the line connectingtheir centers.

TABLE 1 Energy (Method according to an exemplary Data Set embodiment)Energy (J-splines) Spiral (FIG. 1) 12387 12281 Four balls (FIG. 4) 67012386 Mixed overlap (FIG. 5) 1400 3008 M shape (FIG. 6) 1287 34200 Sinewave (FIG. 7) 11447 11441

Table 1 presents the results of the comparison. For each skin, wecomputed the energy J as described in Equation 22. While the result forthe spiral and sine wave was slightly lower for the J-spline method, thedifference is not significant (i.e., the difference is less than 1%) inboth cases and the surfaces have an identical appearance. The methodaccording to an exemplary embodiment of the present inventiondemonstrates significant improvement however for the other surfaces. Theprimary reason for this is that with J-splines, a self-intersectionoccurs that results in high local curvatures and a fatter object ingeneral, as demonstrated in FIG. 9. The method according to an exemplaryembodiment of the present invention will penalize such localself-intersections and deform the spline surface so that the individualsplines are smooth upon convergence.

In other words, as can be seen from FIG. 9, the method according to anexemplary embodiment of the present invention produces smoother resultsthan J-splines. One the left, FIGS. 9( a), (c) and (e), is the J-splineskin for the example shown in FIG. 6. There are two intersections: oneis shown in a zoomed view of FIG. 9( c), with the circle of contactremoved for clarity. In FIG. 9( a), we show the wireframe model of thesplines, and in FIG. 9( c) a texture-mapped view. One the right, FIGS.9( b), (d) and (f) is the result of the method according to an exemplaryembodiment of the present invention, which does not haveself-intersection.

We note that our gradient descent approach guarantees a locally optimalsolution; the particular solution depends on the convexity of the energyfunctional as well as the initial condition. In the examples shown inthis disclosure, the initial skins are chosen to be far from the finalsolution to demonstrate the effect and robustness of the differentialequations. FIG. 10( a) and a zoom in view FIG. 10( b), show ademonstration of ball skinning according to an exemplary embodiment ofthe present invention with a very poor initialization that has a strongself-intersection. Despite the undesirable initialization, the algorithmis able to untangle the self-intersection and produce a smoothinterpolation of the balls, shown in FIG. 10( c). In FIG. 10( d) and azoom in view FIG. 10( e), we show a demonstration of a severeself-intersection where the surface has completely folded in on itself.This initialization is not in the basin of attraction of the desiredsolution, so the skin upon convergence, shown in FIG. 10( f), is not thedesired solution. In practice, it is typically easy to determine a goodinitialization by choosing an initialization for each ball such that thenormal of the intersection plane points along the vector that connectsadjacent ball centroids.

Note that the skin generated by a method according to an exemplaryembodiment of the present invention may pass through a ball (shown inFIG. 11( b)), since it is only constrained to be tangent to the ball atthe circle of intersection. For points not in the circle's plane, theskin may be larger or smaller than the ball. Thus, the skin does notprovide an exact envelope of the balls, but rather, an approximatingenvelope. In an application of modeling blood vessels, this is anacceptable solution since the ball itself is a geometric proxy of thelocal vessel geometry.

In this disclosure, we presented a method for optimally skinning anordered set of 3D balls. Our formulation of the problem requires thatthe skin be modeled by a circle of contact with each ball and that theskin be tangent to the ball along this circle. We have presented novelderivations resulting in differential equations that minimize a convexcombination of the surface area and mean curvature of a third orderpolynomial spline surface subject to these constraints. Starting with aninitial skin, we evolve the skin's parameters until convergence.Experimental results demonstrate the viability of this method.

Derivatives of the coefficients of the first fundamental form are nowprovided.

The derivatives for the j the ball's (j ε[i−1,i]) coefficients of thefirst fundamental form with respect to w^(k) εw^(i)=[θ^(i),φ^(i),d^(i)]^(T) are given by

$\frac{\partial E^{j}}{\partial w^{k}} = {2\; {S_{u}^{j} \cdot \frac{\partial S_{u}^{j}}{\partial w^{k}}}}$$\frac{\partial F^{j}}{\partial w^{k}} = {{S_{u}^{j} \cdot \frac{\partial S_{v}^{j}}{\partial w^{k}}} + {S_{v}^{j} \cdot \frac{\partial S_{u}^{j}}{\partial w^{k}}}}$$\frac{\partial G^{j}}{\partial w^{k}} = {2\; {S_{v}^{j} \cdot \frac{\partial S_{v}^{j}}{\partial w^{k}}}}$where S_(u)^(j) = A_(u)^(j)v³ + B_(u)^(j)v² + C_(u)^(j)v + D_(u)^(j)$\frac{\partial S_{u}^{j}}{\partial w^{k}} = {{\frac{\partial A_{u}^{j}}{\partial w^{k}}v^{3}} + {\frac{\partial B_{u}^{j}}{\partial w^{k}}v^{2}} + {\frac{\partial C_{u}^{j}}{\partial w^{k}}v} + \frac{\partial D_{u}^{j}}{\partial w^{k}}}$S_(v)^(j) = 3 A^(j)v² + 2 B^(j)v + C^(j)$\frac{\partial S_{v}^{j}}{\partial w^{k}} = {{3\frac{\partial A^{j}}{\partial w^{k}}v^{2}} + {2\frac{\partial B^{j}}{\partial w^{k}}v} + \frac{\partial C^{j}}{\partial w^{k}}}$and $\frac{\partial A^{j}}{\partial w^{k}} = \left\{ {{\begin{matrix}{{{2\frac{\partial p^{j}}{\partial w^{k}}} + {t^{j}\frac{\partial N^{j}}{\partial w^{k}}}},} & {j = i} \\{{{{- 2}\frac{\partial p^{j + 1}}{\partial w^{k}}} + {t^{j + 1}\frac{\partial N^{j + 1}}{\partial w^{k}}}},} & {j = {i - 1}}\end{matrix}\frac{\partial B^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{{{{- 3}\frac{\partial p^{j}}{\partial w^{k}}} - {2\; t^{j}\frac{\partial N^{j}}{\partial w^{k}}}},} & {j = i} \\{{{3\frac{\partial p^{j + 1}}{\partial w^{k}}} - {t^{j + 1}\frac{\partial N^{j + 1}}{\partial w^{k}}}},} & {j = {i - 1}}\end{matrix}\frac{\partial C^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{{t^{j}\frac{\partial N^{j}}{\partial w^{k}}},} & {j = i} \\{0,} & {j = {i - 1}}\end{matrix}\frac{\partial D^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{\frac{\partial p^{j}}{\partial w^{k}},} & {j = i} \\{0,} & {j = {i - 1}}\end{matrix}A_{u}^{j}} = {{{{- 2}\; p_{u}^{j + 1}} + {2\; p_{u}^{j}B_{u}^{j}}} = {{{3\; p_{u}^{j + 1}} - {3\; p_{u}^{j}C_{u}^{j}}} = {{0D_{u}^{j}} = {{p_{u}^{j}\frac{\partial A_{u}^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{{2\frac{\partial p_{u}^{j}}{\partial w^{k}}},} & {j = i} \\{{{- 2}\frac{\partial p_{u}^{j + 1}}{\partial w^{k}}},} & {j = {i - 1}}\end{matrix}\frac{\partial B_{u}^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{{{- 3}\frac{\partial p_{u}^{j}}{\partial w^{k}}},} & {j = i} \\{{3\frac{\partial p_{u}^{j + 1}}{\partial w^{k}}},} & {j = {i - 1}}\end{matrix}\frac{\partial C_{u}^{j}}{\partial w^{k}}} = {{0\frac{\partial D_{u}^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{\frac{\partial p_{u}^{j}}{\partial w^{k}},} & {j = i} \\{0,} & {j = {i - 1}}\end{matrix}\frac{\partial p^{i}}{\partial\theta^{i}}} = {{a^{i}\frac{\partial R^{i}}{\partial\theta^{i}}\left( {{x\; \cos \; u} + {y\; \sin \; u}} \right)\frac{\partial p^{i}}{\partial\varphi^{i}}} = {{a^{i}\frac{\partial R^{i}}{\partial\varphi^{i}}\left( {{x\; \cos \; u} + {y\; \sin \; u}} \right)p_{u}^{i}} = {{a^{i}{R^{i}\left( {{{- x}\; \sin \; u} + {y\; \cos \; u}} \right)}\frac{\partial p_{u}^{i}}{\partial\theta^{i}}} = {{a^{i}\frac{\partial R^{i}}{\partial\theta^{i}}\left( {{{- x}\; \sin \; u} + {y\; \cos \; u}} \right)\frac{\partial p_{u}^{i}}{\partial\varphi^{i}}} = {{a^{i}\frac{\partial R^{i}}{\partial\varphi^{i}}\left( {{{- x}\; \sin \; u} + {y\; \cos \; u}} \right){and}\frac{\partial R^{i}}{\partial\theta^{i}}} = {{\begin{bmatrix}{{- \sin}\; \theta^{i}\cos \; \varphi^{i}} & {{- \cos}\; \theta^{i}} & {{- \sin}\; \theta^{i}\sin \; \varphi^{i}} \\{\cos \; \theta^{i}\cos \; \varphi^{i}} & {{- \sin}\; \theta^{i}} & {\cos \; \theta^{i}\sin \; \varphi^{i}} \\0 & 0 & 0\end{bmatrix}\frac{\partial R^{i}}{\partial\varphi^{i}}} = {{\begin{bmatrix}{{- \cos}\; \theta^{i}\sin \; \varphi^{i}} & 0 & {\cos \; \theta^{i}\cos \; \varphi^{i}} \\{{- \sin}\; \theta^{i}\sin \; \varphi^{i}} & 0 & {\sin \; \theta^{i}\cos \; \varphi^{i}} \\{{- \cos}\; \varphi^{i}} & 0 & {{- \sin}\; \varphi^{j}}\end{bmatrix}{and}\frac{\partial N^{i}}{\partial\theta^{i}}} = {{\left\lbrack {{{- \sin}\; \theta^{i}\sin \; \varphi^{i}},{\cos \; \theta^{i}\sin \; \varphi^{i}},0} \right\rbrack^{T}\frac{\partial N^{i}}{\partial\varphi^{i}}} = \left\lbrack {{\cos \; \theta^{i}\cos \; \varphi^{i}},{\sin \; \theta^{i}\cos \; \varphi^{i}},{{- \sin}\; \varphi^{i}}} \right\rbrack^{T}}}}}}}}}} \right.}} \right.} \right.}}}}} \right.} \right.} \right.} \right.$

Derivatives of the coefficients of the second fundamental form are nowprovided.

The derivatives for the j the ball's (j ε[i−1,i]) coefficients of thesecond fundamental form with respect to w^(k) εw^(i)=[θ^(i),φ^(i),d^(i)]^(T) are given by

$\frac{\partial e^{j}}{\partial w^{k}} = {{\frac{\partial M^{j}}{\partial w^{k}} \cdot S_{uu}^{j}} + {M^{j} \cdot \frac{\partial S_{uu}^{j}}{\partial w^{k}}}}$$\frac{\partial f^{j}}{\partial w^{k}} = {{\frac{\partial M^{j}}{\partial w^{k}} \cdot S_{uv}^{j}} + {M^{j} \cdot \frac{\partial S_{uv}^{j}}{\partial w^{k}}}}$$\frac{\partial g^{j}}{\partial w^{k}} = {{\frac{\partial M^{j}}{\partial w^{k}} \cdot S_{vv}^{j}} + {M^{j} \cdot \frac{\partial S_{vv}^{j}}{\partial w^{k}}}}$where$M^{j} = \frac{S_{u}^{j} \times S_{v}^{j}}{{S_{u}^{j} \times S_{v}^{j}}}$$\frac{\partial M^{j}}{\partial w^{k}} = {\frac{\frac{\partial S_{u}^{j}}{\partial w^{k}} \times S_{v}^{j}}{{S_{u}^{j} \times S_{v}^{j}}} + \frac{S_{u}^{j} \times \frac{\partial S_{v}^{j}}{\partial w^{k}}}{{S_{u}^{j} \times S_{v}^{j}}} - \frac{\left\lbrack {\left( {S_{u}^{j} \times S_{v}^{j}} \right) \cdot \left( {{\frac{\partial S_{u}^{j}}{\partial w^{k}} \times S_{v}^{j}} + {S_{u}^{j} \times \frac{\partial S_{v}^{j}}{\partial w^{k}}}} \right)} \right\rbrack S_{u}^{j} \times S_{v}^{j}}{{{S_{u}^{j} \times S_{v}^{j}}}^{3}}}$

where

$S_{u}^{j},S_{v}^{j},\frac{\partial S_{u}^{j}}{\partial w^{k}},{{and}\mspace{14mu} \frac{\partial S_{v}^{j}}{\partial w^{k}}}$

are given in Appendix A, and

S_(uu)^(j) = A_(uu)^(j)v³ + B_(uu)^(j)v² + C_(uu)^(j)v + D_(uu)^(j)S_(uv)^(j) = 3 A_(u)^(j)v² + 2 B_(u)^(j)v + C_(u)^(j)S_(vv)^(j) = 6 A^(j)v + 2 B^(j)$\frac{\partial S_{uu}^{j}}{\partial w^{k}} = {{\frac{\partial A_{uu}^{j}}{\partial w^{k}}v^{3}} + {\frac{\partial B_{uu}^{j}}{\partial w^{k}}v^{2}} + {\frac{\partial C_{uu}^{j}}{\partial w^{k}}v} + \frac{\partial D_{uu}^{j}}{\partial w^{k}}}$$\frac{\partial S_{uv}^{j}}{\partial w^{k}} = {{3\frac{\partial A_{u}^{j}}{\partial w^{k}}v^{2}} + {2\frac{\partial B_{u}^{j}}{\partial w^{k}}v} + \frac{\partial C_{u}^{j}}{\partial w^{k}}}$$\frac{\partial S_{vv}^{j}}{\partial w^{k}} = {{6\frac{\partial A^{j}}{\partial w^{k}}v} + {2\frac{\partial B^{j}}{\partial w^{k}}}}$and$\frac{\partial A_{uu}^{j}}{\partial w^{k}} = \left\{ {{\begin{matrix}{{2\frac{\partial p_{uu}^{j}}{\partial w^{k}}},} & {j = i} \\{{{- 2}\frac{\partial p_{uu}^{j + 1}}{\partial w^{k}}},} & {j = {i - 1}}\end{matrix}\frac{\partial B_{uu}^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{{{- 3}\frac{\partial p_{uu}^{j}}{\partial w^{k}}},} & {j = i} \\{{3\frac{\partial p_{uu}^{j + 1}}{\partial w^{k}}},} & {j = {i - 1}}\end{matrix}\frac{\partial C_{uu}^{j}}{\partial w^{k}}} = {{0\frac{\partial D_{uu}^{j}}{\partial w^{k}}} = \left\{ {{\begin{matrix}{\frac{\partial p_{uu}^{j}}{\partial w^{k}},} & {j = i} \\{0,} & {j = {i - 1}}\end{matrix}{where}p_{uu}^{i}} = {{a^{i}{R^{i}\left( {{{- x}\; \cos \; u} - {y\; \sin \; u}} \right)}\frac{\partial p_{uu}^{i}}{\partial\theta^{i}}} = {{a^{i}\frac{\partial R^{i}}{\partial\theta^{i}}\left( {{{- x}\; \cos \; u} - {y\; \sin \; u}} \right)\frac{\partial p_{uu}^{i}}{\partial\varphi^{i}}} = {a^{i}\frac{\partial R^{i}}{\partial\varphi^{i}}\left( {{{- x}\; \cos \; u} - {y\; \sin \; u}} \right)}}}} \right.}} \right.} \right.$

FIG. 12 is a flowchart illustrating a method of ball skinning accordingto an exemplary embodiment of the present invention. As shown in FIG.12, an initial skin S(u,v) is generated (1210). Differential Equation 13is solved to minimize the initial skin's surface area (1220).Differential Equation 21 is solved to minimize a squared mean curvatureof the initial skin's surface (1230). Steps 1220 and 1230 give anupdated skin. Steps 1220 and 1230 are repeated using each newlygenerated updated skin until a desired skin is realized (1240). It isnoted that the updated skin can be achieved by only solving one of thedifferential equations. It is further noted that when both differentialequations are used to generate the updated skin, Equation 23 is used tocombine the differential equations. It is further noted that the numberof iterations can be determined based on a gradient descent procedurerepresented by Equation 24.

A system in which exemplary embodiments of the present invention may beimplemented will now be described with reference to FIG. 13. As shown inFIG. 13, the system includes a scanner 1305, a computer 1315 and adisplay 1310 connected over a wired or wireless network 1320. Thescanner 1305 may be an image scanner that uses a charge-coupled device(CCD) or Contact Image Sensor (CIS), or it may be a magnetic resonance(MR) or computed tomography (CT) scanner, for example. The computer 1315includes, inter alia, a central processing unit (CPU) 1325, a memory1330 and a ball skinning module 1335 that includes program code forexecuting methods in accordance with exemplary embodiments of thepresent invention. The display 1310 is a computer screen, for example.

It is understood that the present invention may be implemented invarious forms of hardware, software, firmware, special purposeprocessors, or a combination thereof. In one embodiment, the presentinvention may be implemented in software as an application programtangibly embodied on a program storage device (e.g., magnetic floppydisk, RAM, CD ROM, DVD, ROM, and flash memory). The application programmay be uploaded to, and executed by, a machine comprising any suitablearchitecture.

It is also understood that because some of the constituent systemcomponents and method steps depicted in the accompanying figures may beimplemented in software, the actual connections between the systemcomponents (or the process steps) may differ depending on the manner inwhich the present invention is programmed. Given the teachings of thepresent invention provided herein, one of ordinary skill in the art willbe able to contemplate these and similar implementations orconfigurations of the present invention.

It is further understood that the above description is onlyrepresentative of illustrative embodiments. For convenience of thereader, the above description has focused on a representative sample ofpossible embodiments, a sample that is illustrative of the principles ofthe invention. The description has not attempted to exhaustivelyenumerate all possible variations. That alternative embodiments may nothave been presented for a specific portion of the invention, or thatfurther undescribed alternatives may be available for a portion, is notto be considered a disclaimer of those alternate embodiments. Otherapplications and embodiments can be implemented without departing fromthe spirit and scope of the present invention.

It is therefore intended, that the invention not be limited to thespecifically described embodiments, because numerous permutations andcombinations of the above and implementations involving non-inventivesubstitutions for the above can be created, but the invention is to bedefined in accordance with the claims that follow. It can be appreciatedthat many of those undescribed embodiments are within the literal scopeof the following claims, and that others are equivalent.

1. A method of computing a continuous interpolation of a discrete set ofthree-dimensional (3D) balls, comprising: generating an initial skin,wherein the initial skin is a surface comprised of splines and whereinthe splines touch each ball along a circle that is tangent to the ball;solving a first differential equation to minimize the initial skin'ssurface area or solving a second differential equation to minimize asquared mean curvature of the initial skin's surface, wherein the resultof solving the first or second differential equations is an updatedskin; and repeating the steps of solving the first or seconddifferential equations for the updated skin, and then, repeating thesteps of solving the first or second differential equations for eachsubsequently updated skin until a desired skin is realized.
 2. Themethod of claim 1, where a circle on a ball B^(i) is in a plane withnormal N^(i)=[cos θ^(i) sin φ^(i), sin φ^(i) sin φ^(i), cos φ^(i)]^(T)that passes through the ball center c^(i) and intersects the ball B^(i),q(u) is a point on the circle and is defined by q(u)={circumflex over(x)} cos u+ŷ sin u, and a point p^(i)(u) on a circle through which aspline passes is defined by p^(i)(u)=c^(i)+r^(i)R^(i)q(u), where r^(i)is a radius of an ith ball, and R^(i) is a rotation matrix specified byN^(i).
 3. The method of claim 2, wherein the rotation matrix R^(i) is$R^{i} = {{R\left( {\theta^{i},\varphi^{i}} \right)} = {\begin{bmatrix}{\cos \; \theta^{i}\cos \; \varphi^{i}} & {{- \sin}\; \theta^{i}} & {\cos \; \theta^{i}\sin \; \varphi^{i}} \\{\sin \; \theta^{i}\cos \; \varphi^{i}} & {\cos \; \theta^{i}} & {\sin \; \theta^{i}\sin \; \varphi^{i}} \\{{- \sin}\; \varphi^{i}} & 0 & {\cos \; \varphi^{i}}\end{bmatrix}.}}$
 4. The method of claim 2, wherein the initial skin isrepresented by S(u,v), where S(u,v) is a collection of continuoussegments S^(i)(u,v), where an S^(i)(u,v) is specified byA^(i)(u),B^(i)(u),C^(i)(u), D^(i)(u), whereA ^(i)(u)=−2p ^(i+1)(u)+2p ^(i)(u)+t ^(i) N ^(i) +t ^(i+1) N ^(i+1),B^(i)(u)=3p ^(i+1)(u)−3p ^(i)(u)−2t ^(i) N ^(i) −t ^(i+1) N ^(i+1),C ^(i)(u)=t ^(i) N ^(i), and D^(i)(u)=p^(i)(u).
 5. The method of claim4, wherein the first differential equation is represented by$\frac{\partial J^{a}}{\partial w^{k}} = {{\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\; F^{i}\frac{\partial F^{i}}{\partial w^{k}}}}{\left\lbrack {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{v}}}}} + {\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\; F^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}}}{\left\lbrack {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{{v}.}}}}}}$6. The method of claim 5, where J^(a)=∫∫√{square root over (EG−F²)}dudv,w^(k) ε[θ^(i),φ^(i)]∀i, E=S_(u)·S_(u), F=S_(u)·S_(v), and G=S_(v)·S_(v).7. The method of claim 4, wherein the second differential equation isrepresented by$\frac{\partial J^{c}}{\partial w^{k}} = {{\int{\int{2\; {H^{i}\left\lbrack {\frac{{\frac{\partial\varepsilon^{i}}{\partial w^{k}}G^{i}} + {e^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial f^{i}}{\partial w^{k}}F^{i}} - {2\; f^{i}\frac{\partial F^{i}}{\partial w^{k}}} + {\frac{\partial g^{i}}{\partial w^{k}}E^{i}} + {g^{i}\frac{\partial E^{i}}{\partial w^{k}}}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)} - \frac{\left( {{\varepsilon^{i}G^{i}} - {2\; f^{i}F^{i}} + {g^{i}E^{i}}} \right) \cdot \left( {{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial F^{i}}{\partial w^{k}}F^{i}}} \right)}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)^{2}}} \right\rbrack}{u}{v}}}} + {\int{\int{2\; {H^{i - 1}\left\lbrack {\frac{\begin{matrix}{{\frac{\partial\varepsilon^{i - 1}}{\partial w^{k}}G^{i - 1}} + {e^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial f^{i - 1}}{\partial w^{k}}F^{i - 1}} -} \\{{2\; f^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}} + {\frac{\partial g^{i - 1}}{\partial w^{k}}E^{i - 1}} + {g^{i - 1}\frac{\partial E^{i - 1}}{\partial w^{k}}}}\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)} - \frac{\begin{matrix}{\left( {{e^{i - 1}G^{i - 1}} - {2\; f^{i - 1}F^{i - 1}} + {g^{i - 1}E^{i - 1}}} \right) \cdot} \\\left( {{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial F^{i - 1}}{\partial w^{k}}F^{i - 1}}} \right)\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)^{2}}} \right\rbrack}{u}{{v}.}}}}}$8. The method of claim 7, where J^(c)=∫∫H²dudv, w^(k) ε[θ_(i),φ^(i)]∀i,E=S_(u)·S_(u), F=S_(u)·S_(v), G=S_(v)·S_(v), e=M·S_(uu, f=M·S) _(uv),and g=M·S_(vv), where M is the surface normal.
 9. The method of claim 4,further comprising: combining the first and second differentialequations to generate the updated skin, wherein the first and seconddifferential equations are combined by the following equation${\frac{\partial J}{\partial w^{k}} = {{\left( {1 - k} \right)\frac{\partial J^{a}}{\partial w^{k}}} + \frac{\partial J^{c}}{\partial w^{k}}}},$where $\frac{\partial J^{a}}{\partial w^{k}}$ is the first differentialequation and $\frac{\partial J^{c}}{\partial w^{k}}$ is the seconddifferential equation.
 10. The method of claim 9, wherein the firstdifferential equation is represented by$\frac{\partial J^{a}}{\partial w^{k}} = {{\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\; F^{i}\frac{\partial F^{i}}{\partial w^{k}}}}{\left\lbrack {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{v}}}}} + {\frac{1}{2}{\int{\int{\frac{{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\; F^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}}}{\left\lbrack {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right\rbrack^{\frac{1}{2}}}{u}{{v}.}}}}}}$11. The method of claim 10, where J^(a)=∫∫√{square root over(EG−F²)}dudv, w^(k) ε[θ^(i),φ^(i)]∀i, E=S_(u)·S_(u), F=S_(u)·S_(v), andG=S_(v)·S_(v).
 12. The method of claim 9, wherein the seconddifferential equation is represented by$\frac{\partial J^{c}}{\partial w^{k}} = {{\int{\int{2\; {H^{i}\left\lbrack {\frac{{\frac{\partial\varepsilon^{i}}{\partial w^{k}}G^{i}} + {e^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial f^{i}}{\partial w^{k}}F^{i}} - {2\; f^{i}\frac{\partial F^{i}}{\partial w^{k}}} + {\frac{\partial g^{i}}{\partial w^{k}}E^{i}} + {g^{i}\frac{\partial E^{i}}{\partial w^{k}}}}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)} - \frac{\left( {{e^{i}G^{i}} - {2\; f^{i}F^{i}} + {g^{i}E^{i}}} \right) \cdot \left( {{\frac{\partial E^{i}}{\partial w^{k}}G^{i}} + {E^{i}\frac{\partial G^{i}}{\partial w^{k}}} - {2\frac{\partial F^{i}}{\partial w^{k}}F^{i}}} \right)}{2\left( {{E^{i}G^{i}} - \left( F^{i} \right)^{2}} \right)^{2}}} \right\rbrack}{u}{v}}}} + {\int{\int{2\; {H^{i - 1}\left\lbrack {\frac{\begin{matrix}{{\frac{\partial\varepsilon^{i - 1}}{\partial w^{k}}G^{i - 1}} + {e^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial f^{i - 1}}{\partial w^{k}}F^{i - 1}} -} \\{{2\; f^{i - 1}\frac{\partial F^{i - 1}}{\partial w^{k}}} + {\frac{\partial g^{i - 1}}{\partial w^{k}}E^{i - 1}} + {g^{i - 1}\frac{\partial E^{i - 1}}{\partial w^{k}}}}\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)} - \frac{\begin{matrix}{\left( {{e^{i - 1}G^{i - 1}} - {2\; f^{i - 1}F^{i - 1}} + {g^{i - 1}E^{i - 1}}} \right) \cdot} \\\left( {{\frac{\partial E^{i - 1}}{\partial w^{k}}G^{i - 1}} + {E^{i - 1}\frac{\partial G^{i - 1}}{\partial w^{k}}} - {2\frac{\partial F^{i - 1}}{\partial w^{k}}F^{i - 1}}} \right)\end{matrix}}{2\left( {{E^{i - 1}G^{i - 1}} - \left( F^{i - 1} \right)^{2}} \right)^{2}}} \right\rbrack}{u}{{v}.}}}}}$13. The method of claim 12, where J^(c)=∫∫H²dudv, w^(k)ε[θ^(i),φ^(i)]∀i, E=S_(u)·S_(u), F=S_(u)·S_(v), G=S_(v)·S_(v),e=M·S_(uu), f=M·S_(uv), and g=M·S_(vv), where M is the surface normal.14. The method of claim 9, wherein the first and second differentialequations are used in a gradient descent procedure to find the desiredskin.
 15. The method of claim 14, wherein the gradient descent proceduremanipulates parameters w^(i)=[θ^(i),φ^(i)]^(T) of each ball i, wherew^(i)(n+1)=w^(i)(n)−Δt∇J_(w) _(i) _((n).)∀i to find the desired skin.16. A method of modeling a tubular structure, comprising: imaging atubular structure; placing a plurality of balls in the tubularstructure; and finding a skin that smoothly interpolates the balls,wherein finding the skin comprises: generating an initial skin, whereinthe initial skin is a surface comprised of splines and wherein thesplines touch each ball along a circle that is tangent to the ball;solving a first differential equation to minimize the initial skin'ssurface area or solving a second differential equation to minimize asquared mean curvature of the initial skin's surface, wherein the resultof solving the first or second differential equations is an updatedskin; and repeating the steps of solving the first or seconddifferential equations for the updated skin, and then, repeating thesteps of solving the first or second differential equations for eachsubsequently updated skin until the skin that smoothly interpolates theballs is realized, wherein the skin that smoothly interpolates the ballsis a model of the tubular structure.
 17. The method of claim 16, whereinthe tubular structure is an anatomical structure.
 18. The method ofclaim 16, wherein the tubular structure is imaged by a scanner.
 19. Themethod of claim 16, wherein the balls are ordered.
 20. A system ofcomputing a continuous interpolation of a discrete set ofthree-dimensional (3D) balls, comprising: a memory device for storing aprogram; and a processor in communication with the memory device, theprocessor operative with the program to perform a method, the methodcomprising: generating an initial skin, wherein the initial skin is asurface comprised of splines and wherein the splines touch each ballalong a circle that is tangent to the ball; solving a first differentialequation to minimize the initial skin's surface area or solving a seconddifferential equation to minimize a squared mean curvature of theinitial skin's surface, wherein the result of solving the first orsecond differential equations is an updated skin; and repeating thesteps of solving the first or second differential equations for theupdated skin, and then, repeating the steps of solving the first orsecond differential equations for each subsequently updated skin until adesired skin is realized.